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1  Final  Report 


This  report  documents  the  major  project  findings  during  the  duration  of  the  AFOSR  Young  Investigator 
Program  titled  Information  Collection  and  Fusion  for  Space  Situation  Awareness.  The  principal  goals  of 
this  on-going  research  work  are  threefold: 

1 .  To  understand  how  uncertain  input  variables  such  as  geometric  properties  of  RSO  and  atmospheric 
drag  or  solar  radiation  pressure  (SRP)  coefficient  and  random  stochastic  forcing  affect  the  output  of  a 
dynamical  model  for  orbit  propagation. 

2.  To  estimate  orbit  states  together  with  quantitative  measures  of  confidence  associated  with  those  esti¬ 
mates  by  combining  system  prediction  with  observations  from  various  resources. 

3.  To  design  robust  methodology  for  optimal  sensor  management  while  taking  into  account  the  uncer¬ 
tainties  in  the  system  dynamics. 

The  highlights  of  the  achievements  of  this  research  work  are  described  as  follows: 


2  Adaptive  Gaussian  Mixture  Model: 

A  major  part  of  this  research  work  focused  on  developing  an  approach  to  improve  the  statistical  validity 
of  Resident  Space  Objects  (RSO)  orbit  estimates  and  uncertainties  as  well  as  a  method  of  associating  obser¬ 
vations  with  the  correct  RSO  and  classifying  events  in  near  real  time.  Our  approach  involved  studying  an 
adaptive  Gaussian  mixture  model  (AGMM)  solution  to  the  Fokker-Planck-Kolmogorov  Equation  (FPKE) 
for  its  applicability  to  the  RSO  tracking  problem.  The  FPKE  describes  the  time-evolution  of  stochastic 
systems  and  the  Adaptive  Gaussian  Sum  Filter  (AGSF)  is  a  solution  to  the  equation  that  allows  for  non- 
Gaussian  pdfs.  The  AGSF  algorithm  is  designed  to  be  scalable,  relatively  efficient  for  solutions  of  this 
type,  and  able  to  handle  the  nonlinear  effects  which  are  common  in  the  estimation  of  RSO  orbit  states.  In 
addition,  techniques  for  data  association  that  are  compatible  with  the  AGSF  based  on  entropy  theory  were 
examined.  The  AGSF  and  corresponding  observation  association  methods  were  evaluated  using  simulated 
data  to  determine  their  performance  and  feasibility.  The  key  idea  of  the  AGSF  is  to  approximate  the  state 
pdf  by  a  finite  sum  of  Gaussian  density  functions  whose  mean  and  covariance  are  propagated  from  one 
time-step  to  the  next  using  linear  theory.  The  weights  of  the  Gaussian  kernels  are  updated  at  every  time-step 
by  requiring  the  sum  to  satisfy  the  FPKE.  When  properly  formulated,  the  mixture  problem  in  the  AGSF  can 
be  solved  efficiently  and  accurately  using  convex  optimization  solvers,  even  if  the  mixture  model  includes 
many  terms.  This  methodology  effectively  decouples  a  large  uncertainty  propagation  problem  into  many 
small  problems.  As  a  consequence,  the  solution  algorithm  can  be  parallelized  on  most  High  Performance 
Computing  (HPC)  systems.  Finally,  a  Bayesian  framework  can  be  used  on  the  AGSF  structure  to  assimilate 
(noisy)  observational  data  with  model  forecasts. 

There  is  no  doubt  that  number  of  Gaussian  kernels  in  the  mixture  model  approximation  plays  an  impor¬ 
tant  role  in  its  accuracy  and  computational  complexity.  Adaptive  Gaussian  Sum  Mixture  (AGMM)  method 
to  solve  the  Fokker-Planck-Kolmogorov  Equation  (FPKE)  has  been  modified  to  automatically  select  number 
of  Gaussian  kernels  based  upon  FPKE  error  feedback.  Critical  Gaussian  kernels  are  identified  which  needs 
to  be  split  and  merged  for  better  mixture  approximation.  Furthermore,  sparse  approximation  tools  norm 
minimization)  have  been  used  to  obtain  a  compact  representation  of  the  mixture  model. 
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Figure  [I]  shows  the  contours  of  the  RSO  orbital  state  pdf  at  different  time.  The  RSO  orbit  corresponds  to 
low  Earth  orbit  at  an  altitude  of  700km.  Atmospheric  drag  and  first  harmonic  of  non-sphercal  gravity  Initial 
pdf  was  approximated  using  one  Gaussian  component  having  mean  /x0  and  covariance  So. 

ju0  =  [7000  0  0  0  -  1.0374  7.4771]  (1) 

S0  =  diag  (0.01,  0.01,  0.01,  0.000001,  0.000001, 0.000001)  (2) 

The  initial  mean  and  covariance  were  propagated  for  6  hours  using  the  equations  of  unscented  Kalman 
filter.  Number  of  Gaussian  kernels  in  mixture  approximation  were  automatically  selected  by  making  use  of 
adaptive  split  and  merge  strategy  in  conjunction  with  Kullbeck-Leibler  measure.  125000  monte  carlo  runs 
were  used  as  truth  to  verify  the  developed  approach.  Figure  [I]  shows  the  contours  of  the  pdf  at  different 
time.  One  can  see  that  the  initial  Gaussian  pdf  characterised  by  circular  contour  evolves  into  non-Gaussian 
pdf  over  time.  The  contours  of  the  pdf  match  with  the  contours  of  the  monte  carlo  run.  One  can  also  see  the 
change  in  the  number  of  components  at  different  time  as  a  result  of  split  and  merge  technique. 

The  findings  of  this  work  are  disseminated  through  following  publications: 

•  K.  Vishwajeet,  P.  Singla  and  M.  Jah,  “Nonlinear  Uncertainty  Propagation  for  Perturbed  Two-Body  Or¬ 
bits,”  AIAA  Journal  of  Guidance,  Control  and  Dynamics ,  Accepted,  January  2014,  DOI:  10.25 14/1  .G000472 , 

•  K.  Vishwajeet  and  P.  Singla,  “Adaptive  Split  and  Merge  Technique  for  Gaussian  Mixture  Models  to 
Solve  Kolmogorov  Equation,”  2014  American  Control  Conference,  Portland,  OR,  June  4-6,  2014. 

•  K.  Vishwajeet  and  P.  Singla,  “Adaptive  Split  and  Merge  Algorithm  for  Gaussian  Mixture  Models,” 

2013  Astrodynamics  Specialist  Conference,  Hilton  Head,  North  Carolina. 

•  K.  Vishwajeet  and  P.  Singla,  “Sparse  Approximation  Based  Gaussian  Mixture  Model  Approach  for 
Uncertainty  Propagation  for  Nonlinear  Systems,”  2013  American  Control  Conference,  Washington 
D.C. 

•  K.  Vishwajeet  “Adaptive  Gaussian  Mixture  Model  for  Uncertainty  Propagation  Through  Perturbed 
Two-Body  Model,”  M.S.  thesis,  Department  of  Mechanical  &  Aerospace  Engineering,  University  at 
Buffalo,  Buffalo,  NY,  August  2013.  (Best  M.S.  Thesis  Award  from  NAGS). 

•  G.  Terejanu,  P.  Singla,  T.  Singh  and  P.  Scott,  “Adaptive  Gaussian  Sum  Filter  for  Nonlinear  Bayesian 
Estimation,”  IEEE  Transactions  on  Automatic  Control,  Vol.  56,  Issue  9,  pp.  2151-2156,  Sep.  2011, 

DOI:  10. 1109/TAC. 201 1.2141550, 

3  Conjugate  Unscented  Transform: 

The  problem  of  uncertainty  characterization  in  nonlinear  systems  subject  to  stochastic  excitation  and 
uncertain  initial  conditions  is  of  central  interest  to  various  domains  of  science  and  engineering.  One  may  be 
interested  in  predicting  the  probability  of  collision  of  an  asteroid  with  Earth,  control  of  movement  and  plan¬ 
ning  of  actions  of  autonomous  systems,  diffusion  of  toxic  materials  through  atmosphere,  the  optimization 
of  financial  investment  policies,  or  simply  the  computation  of  the  prediction  step  in  the  implementation  of  a 
Bayes’  filter.  All  these  applications  require  the  study  of  the  relevant  stochastic  system  and  involve  comput¬ 
ing  multi-dimensional  expected  value  integrals  with  respect  to  an  appropriate  probability  density  function 
(pdf).  Analytical  expressions  for  these  multi-dimension  integrals  in  general  exist  for  linear  systems.  For 
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Figure  1:  Contours  corresponding  to  prior  pdf  in  cartesian  coordinates  with  drag  using  MC(cots)  & 
AGSF(thick  line)  at  different  time(’*’  indicates  the  centres  of  components) 


3 


example,  the  well  celebrated  Kalman  filter  provides  the  analytical  expressions  for  the  mean  and  covariance 
of  the  states  of  the  linear  system  subject  to  Gaussian  white  noise  and  Gaussian  initial  condition  errors.  How¬ 
ever  in  most  other  cases  one  often  does  not  have  direct  analytical  solution  for  these  integrals  and  has  to 
approximate  integral  values  by  making  use  of  computational  methods. 

Several  computation  techniques  exist  in  the  literature  to  approximate  the  expectation  integral  with  respect 
to  a  Gaussian  pdf,  the  most  popular  being  Monte  Carlo  (MC)  methods,  Gauss-Hermite  (GH)  Quadrature 
Rules,  Sparse  Grid  quadratures  such  as  Gauss-Hermite  Smolyak  (GHS)  quadratures,  Unscented  Transfor¬ 
mation  (UT)  and  Cubature  methods.  All  these  methods  involve  an  approximation  of  the  expectation  integral 
as  a  weighted  sum  of  integrand  values  at  specified  points  within  the  domain  of  integration.  These  methods 
basically  differ  from  each  other  in  the  generation  of  these  specific  points.  For  example,  MC  methods  involve 
random  samples  from  the  specified  pdf  while  the  Gaussian  quadrature  scheme  involves  deterministic  points 
carefully  chosen  to  reproduce  exactly  the  integrals  for  polynomials.  Quadrature  rules  in  higher  dimensions 
are  usually  referred  to  as  cubature  rules.  A  cubature  rule  is  said  to  be  exact  to  degree  d,  if  it  can  only  in¬ 
tegrate  all  polynomials  with  degree  <  d.  For  ID  (1 -Dimensional)  integrals,  one  needs  n  quadrature  points 
according  to  the  Gaussian  quadrature  scheme  to  exactly  reproduce  the  expectation  integrals  of  polynomi¬ 
als  with  degree  2n  —  1  or  less.  However,  in  generic  ND,  one  needs  to  take  the  tensor  product  of  n-lB 
quadrature  points  and  hence  would  yield  a  total  of  nN  quadrature  points.  This  cubature  rule  is  referred  to 
as  Gaussian  Quadrature  Product  rule.  Even  for  a  moderate-dimension  system  involving  6  random  variables, 
the  number  of  points  required  to  evaluate  the  expectation  integral  with  only  5  points  along  each  direction 
is  56  =  15,  625.  This  is  a  non-trivial  number  of  points  that  might  make  the  calculation  of  the  integral 
computationally  expensive,  especially  when  the  evaluation  of  function  at  each  cubature  point  itself  can  be 
an  expensive  procedure,  e.g.,  one  may  need  to  solve  a  system  of  partial  differential  equations  to  compute 
the  function  of  interest.  The  sparse  grid  quadratures  or  Smolyak  quadrature  schemes  in  particular  take  the 
sparse  product  of  ID  quadrature  rules  and  thus  have  fewer  points  than  the  equivalent  Gaussian  quadrature 
rules  at  the  cost  of  introducing  negative  weights.  But  fortunately  the  Gaussian  quadrature  rule  is  not  minimal 
for  N  >  2  and  there  exists  cubature  rules  with  reduced  number  of  points.  This  forms  the  motivation  for  the 
work  presented  in  this  paper. 

In  the  perspective  of  nonlinear  filtering  where  the  integrals  with  Gaussian  probability  density  functions 
frequently  arise,  cubature  methods  such  as  the  UKF  with  reduced  number  of  points  have  tremendous  poten¬ 
tial  especially  in  an  online  scenario.  The  highly  acclaimed  Unscented  Transform  (UT)  with  only  2N  +  1 
points  is  a  degree  3  cubature  method,  where  these  cubature  points  were  called  sigma  points.  Similar  to  the 
UT  method,  a  more  recent  development  is  the  Cubature  Kalman  filter  (CKF)  which  is  again  exact  to  degree 
3  but  uses  only  2 N  points.  It  can  be  noticed  that  the  sigma/cubature  points  of  2N  +  1  UT  and  2N  CKF 
can  be  observed  as  special  cases  of  a  more  general  structure  of  cubature  points  previously  presented  in  for 
symmetrical  probability  density  functions. 

In  this  research  work,  a  new  methodology  is  developed  to  efficiently  evaluate  expectation  integrals  in 
general  TV-dimensional  space  by  satisfying  higher  order  moment  constraint  equations.  New  sets  of  cuba- 
ture/quadrature  points  are  defined  to  satisfy  moment  equations  up  to  the  tenth  order.  The  developed  method¬ 
ology  can  be  used  as  an  efficient  alternative  to  Gaussian  quadrature  rule  with  significantly  reduced  number 
of  function  evaluations  while  maintaining  accuracy.  The  proposed  methodology  extends  the  main  idea  of 
the  conventional  unscented  transformation  method  to  construct  a  fully  symmetric  reduced  sigma/cubature 
point  set  with  all  positive  weights  that  sum  up  to  one  and  are  equivalent  to  the  Gaussian  quadrature  product 
rule  of  same  order.  In  a  numerical  context,  equivalent  to  same  order  implies  that  for  a  polynomial  of  order 
2m  —  1  in  N-dimensions,  both  the  new  reduced  sigma  points  from  the  proposed  method  known  as  Conjugate 
Unscented  Transform  method  (CUT)  and  the  mN  quadrature  points  from  the  Gaussian  quadrature  product 
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rule  result  in  same  order  of  relative  percentage  error.  The  main  idea  of  the  CUT  approach  is  to  judiciously 
select  specific  structures  for  sigma  points  rather  than  taking  tensor  product  of  1-D  points  like  in  the  Gauss 
quadrature  scheme.  This  allows  us  to  compute  multi-dimension  expectation  integrals  with  the  same  orders 
of  magnitude  in  terms  of  error  but  achieves  this  with  a  far  less  number  of  points. 


(c)  Number  of  points  required  to  achieve  0.5% 
Rel.  error 


Figure  2:  GH  vs.  CUT:  Simulation  Results  corresponding  to  /(X)  =  (Vl  +  XTX)  3 

Fig.  |2(a)|  shows  the  relative  error  with  respect  to  the  true  value  of  the  following  expectation  integral  as  a 
function  of  GH  quadrature  points  along  each  dimension  while  varying  the  dimension  of  input  space  from  2 
to  6. 


£[/(X)]  = 


+  XTX)“,  X  e  RN  xOatxiO.1  INxNchi 


(3) 


This  is  a  benchmark  problem  and  it  had  been  discussed  that  computing  the  expectation  of  /(X)  for  negative 
values  for  a  is  a  challenging  task  since  a  <  0  leads  to  a  delta- sequence  functions.  Taking,  a  =  —  3  and  the 
identity  covariance  of  the  zero  mean  Gaussian  pdf  is  scaled  down  by  0.1  to  make  the  integral  value  converge 
with  reasonable  number  of  points.  Similarly,  Fig.  |2(b)  shows  the  relative  error  for  dimensions  2  to  6  using 
the  CUT  methods.  Furthermore,  Fig.  2(c)|  shows  the  minimal  number  of  points  required  for  each  method  to 
achieve  at  least  0.5%  relative  error.  It  is  clear  that  the  proposed  CUT8  method  yields  the  accuracy  of  0.5% 
as  the  GH4  quadrature  rule  but  with  much  smaller  fraction  of  the  number  of  points. 

One  could  take  advantage  of  these  cubature  rules  in  the  evaluation  of  the  integrals  while  maintaining 
higher  accuracy  and  thus  feasible  on-line  execution  of  a  filter.  To  illustrate  the  effectiveness  of  the  proposed 
approach,  let  us  consider  a  typical  air  traffic  control  scenario  as  shown  in  Figure  |3(a)|  The  kinematics  of 
the  turning  motion  is  modeled  by  the  set  of  nonlinear  equations  called  ‘CT-  Coordinated  Turn.  The  CT 
model  is  characterized  by  constant  speed  and  constant  turn  rate.  The  turn  rate  Q  is  usually  unknown  and 
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is  hence  appended  to  the  state  vector  making  the  model  nonlinear.  The  system  dynamic  equations  for  CT- 
model,  where  the  state  vector  x  =  [£  £  77  77  fl]T  and  sensor  model  that  consists  of  a  radar  located  at 

the  origin  and  measures  the  range  and  bearing  are: 
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The  process  noise  Vk  and  measurement  noise  uok  are  independent  zero  mean  gaussian  noise  processes  with 
following  covariance  matrices: 
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The  parameters  used  in  the  simulation  are  given  as: 


T  =  5s 


L\  =  0.16  ar  =  100  m 


Tf  =  495s  L2  0.01  <jq  1  degree 

condition  uncertainty  is: 

xo  =  [25000  771,  —  120  7Ti/s,  10000  772,  0  772/s,  0.000001  rad/s]T 
Pq/ 0  =  diag([10002  7722,  100  72i2/s2,  10002  7212,  100  tt22/s2,  (17t/180)2  rad2/s2]) 

The  filters  used  in  the  simulations  are  listed  below: 


The  initial 
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Cubature  Kalman  Filter  (CKF)-  with  10  points 

Unscented  Kalman  Filter  (UKF)  -  with  1 1  points 

CUT4-  with  42  points 

CUT6-  with  83  points 

CUT8-  with  355  points 

GH6  -  with  7776  points 

Particle  Filter  (PF)-  with  5000  sample  points 

- True  Path 

•  Measurements 
■  Sensor 


-1.5  -1  -0.5  0  0.5 


1.5  2  2.5  3 


(a)  Airplane  Trajectory  (b)  Filter  Estimates 

Figure  3:  Simulation  Results  for  Air-Traffic  Scenario 


A  simulation  of  all  these  filters  are  shown  in  Figure  |3(b)|  The  trajectories  of  each  filter  shown  in  the 
Figure [3(b)] is  the  average  over  200  runs. 
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As  the  dimension  of  the  process  model  is  5,  to  avoid  a  negative  weight  in  the  UT  a  kappa  value  of 
ft  =  1  is  considered.  Each  run  consists  of  randomly  selecting  a  point  from  initial  distribution  as  the  initial 
condition  for  the  filters,  generating  the  measurements  from  the  true  trajectory  by  the  sensor  model  ([4])  and 
random  noise  with  covariance  R.  All  the  filters  are  initiated  with  the  same  initial  point  and  use  the  same  set 


of  measurements  for  a  particular  run.  A  batch  of  200  runs  are  performed  for  each  filter  considered.  The  root 
mean  square  error  RMSE  in  position, velocity  and  turn  rate  for  the  batch  of  runs  is  calculated  for  each  filters 
as  follows 

RMSEpos(k)  =  . 

1  200 

200  “  &(*0)2  +  fo(fc)  -  Vi(k))2) 

i=l 

(8) 

RMSEvel(k)  =  . 

200 

200  “  &(*0)2  +  -  Vi(k))2) 

i=l 

(9) 

RMSEQ{k)  =  . 

1  200 

i=l 

(10) 

The  RMSE  in  position  and  velocity  is  calculated  individually  for  each  filter  estimate  with  respect  to  the  true 
trajectory.  The  2-norm  in  the  RMSE  is  calculated  as: 


\\RMSEpos\\2 


N 


nTf 

RM  S  Epos(k)2 , 

k= 1 


\\RMSEvel  ||  2  = 


N 


Y  RMSE„d(kf 

k= 1 


(11) 


where  nrf  is  the  number  time  steps  in  the  simulation.  Table  [I]  gives  a  quick  summary  of  the  2-norms  in  the 
RMSE  errors.  Illustrative  figures  are  also  shown  to  compare  the  various  filters.  The  RMSEpos  is  shown 
in  Figures  [4(a)]-  |4(d)|  For  the  sake  of  readability,  Figure  |4(a)|  shows  the  performance  of  all  the  filters  while 
the  rest  of  the  Figures  |4(b)||4(d)|  show  the  same  results  but  with  different  combinations  of  filters  to  avoid  the 
overlap  of  the  results  of  different  filters.  In  the  figures  ‘PF-rnu’  denotes  the  mean  estimate  of  the  particle 
filter  while  4PF-mo’  represents  the  mode  estimate  of  the  particle  filter.  One  can  see  that  the  estimates  of  the 
Particle  filter  mean  and  mode  are  very  close.  One  can  see  the  gradual  improvement  as  the  order  of  the  filters 
are  increased.  CKF  and  UKF  being  3rd  order  rules  have  the  maximum  error.  While  the  results  of  CUT8  and 
GH6  are  very  much  comparable  to  the  particle  filter. 


Table  1 :  Comparison  of  2-norms  of  RMSE  in  position,  velocity  and  turn  rate 


\\RMSE\\2  in 

PF  —  mean 

PF  —  mode 

CKF 

UKF 

CUT  4 

CUTS 

CUTS 

GH  6 

Position 

2655.80 

2843.22 

23028.35 

16414.90 

4251.00 

3151.78 

3044.16 

2982.91 

Velocity 

544.38 

810.77 

374946.32 

317452.34 

119231.01 

58900.37 

776.35 

895.45 

Cl 

0.8842 

1.6476 

63.4233 

53.7346 

26.4342 

19.2502 

2.0175 

2.2820 

Finally,  this  newly  developed  method  is  tested  for  the  orbit  uncertainty  characterization  for  low  Earth 
object.  For  simulation  purposes,  we  consider  a  satellite  in  Fow  Earth  Orbit  (FEO)  at  an  altitude  of  622  km. 
The  initial  state  probability  density  function  (pdf)  of  the  system  is  assumed  to  be  Gaussian  with  the  following 
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Figure  4:  Comparison  of  RMSEpos  with  Particle  Filter  Results 
mean  (in  km)  and  covariance  (in  km2): 
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The  covariance  matrix  reflects  a  standard  deviation  of  100m  in  the  initial  estimate  of  the  three  position 
coordinates  and  a  standard  deviation  of  1  m/s  in  the  estimate  of  velocity.  The  time  period  of  the  satellite  is 
0.8104  hours.  Exponential  atmospheric  density  model  was  considered  to  evaluate  atmospheric  drag  term. 
Furthermore,  first  harmonic  of  non-spheric  gravity  field,  J2  was  also  considered.  To  assess  the  accuracy 
of  CUT,  we  consider  100,000  Monte  Carlo  (MC)  runs  to  be  ground  truth  and  compare  first  four  moments 
computed  through  CUT  to  with  those  computed  through  MC  runs.  There  are  3,  6,  10  and  15  distinct  first 
order,  second  order,  third  order  and  fourth  order  moments.  At  each  time  time  t ,  the  2-norm  in  the  error  of 
each  moment  is  taken  with  respect  to  the  corresponding  moment  calculated  by  100,000  MC  runs.  These  are 
then  denoted  as  M±(t),  M2(t),  M^(t)  and  M^t).  The  values  are  further  scaled  by  the  radius  of  the  Earth 
and  the  2-norm  over  time  is  taken.  The  results  are  shown  in  Table  [2]  From  these  results,  it  is  clear  that  CUT 
methodology  capture  the  moments  of  non-Gaussian  pdf  with  a  order  of  magnitude  less  points.  This  really 
demonstrate  the  efficacy  of  the  newly  developed  CUT  approach  in  capturing  non-Gaussian  behavior  with 
very  few  model  runs. 

Furthermore,  the  principle  of  maximum  entropy  is  used  to  approximate  orbit  state  pdf  from  moments 
computed  through  CUT  points.  Figure [5]  shows  the  marginalized  pdfs  in  (x,  y )  and  (x,  z)  respectively  using 
the  MC  sample  points  in  3D.  It  can  be  seen  that  by  using  only  moments  upto  order  2,  the  resultant  pdf  is 
Gaussian  that  tends  to  neglect  the  tail  region  of  the  pdf.  As  higher  order  moments  are  considered,  the  tail 
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Table  2:  Error  in  moments  (x,y,z)  with  respect  to  100,000  MC  runs 


Moment 

MC 

UT 

CUT4 

CUT6 

CUT8 

\\Mi{t)/Re\\2 

0.0015 

0.0003 

0.0003 

0.0003 

0.0003 

(||M2(i)/.Re||2)1/2 

0.4802 

0.1635 

0.1650 

0.1650 

0.1650 

(||M3(i)/i?e||2)1/2 

0.7034 

2.0967 

1.3654 

1.3654 

1.3654 

(||M4(t)/i?e||2)1/2 

7.7290 

16.4402 

3.8991 

3.8833 

3.8833 

No.  of  Points 

10000 

13 

76 

137 

745 

region  is  better  captured.  The  marginalized  pdfs  using  CUT4  points  are  shown  in  Figure[6j  The  marginalized 
pdfs  using  CUT6  points  are  shown  in  Figure|7J  The  CUT6  points  show  improvement  compared  to  the  CUT4 
points  and  very  much  the  same  as  the  one  obtained  by  MC  runs.  This  once  again  reiterate  the  significance 
of  newly  developed  CUT  approach  in  capturing  non-Gaussian  behavior. 

Furthermore,  a  test  case  is  considered  to  show  the  efficacy  of  the  proposed  approach  in  approximating 
the  probability  of  collision  between  two  satellites.  The  initial  conditions  are  chosen  such  that  at  the  point  of 
closest  approach,  the  pdf  of  one  satellite  is  significantly  non  Gaussian.  The  satellite  initial  condition  state 


vectors  of  the  form  [x  km ,  y  km ,  z  fcm,  x  km/ s,  y  km/ s,  i  km/ s]  are  taken  as 

»i  =  [7000,  0,  0, 1.0374,  -1.0374,  7.4771]  (12) 

/i2  =  [6729.4309,  -318.1595,  -1381.8922,  2.2649,  -1.1893,  7.1694]  (13) 

Pi  =  [0.01,  0.01,  0.01,  le(“9),  le(“9),  le(“9)]  (14) 

P2  =  [0.01,  0.01,  0.01,  le("8\  le(“8),  le(“8)]  (15) 


The  first  satellite  is  propagated  for  1,  72, 860  seconds  while  the  second  satellite  is  propagated  for  21,  600 
seconds.  The  satellite  that  is  propagated  for  over  two  days  has  significant  non-gaussian  nature.  For  example, 
Figure  [8]  shows  a  snapshot  of  100,000  Monte  Carlo  runs  when  the  random  point  clouds  of  both  satellites 
intersect.  As  the  time  of  closest  approach  is  itself  a  random  variable,  one  does  not  know  the  exact  time 
of  maximum  probability  of  conjunction.  Hence,  the  conjunction  probability  is  evaluated  for  ±20  second 
time  interval  from  the  time  of  closest  approach  calculated  from  the  nominal  (mean)  trajectories  of  the  two 
satellites. 

Figure [9] shows  the  probability  of  collision  computed  through  various  methods.  2.3  million  Monte  Carlo 
(MC)  samples  are  considered  to  provide  ground  truth.  The  full  non  Gaussian  conjunction  probability  com¬ 
puted  via  the  principle  of  maximum  entropy  closely  agrees  with  the  MC  probability.  Further,  the  Gaussian 
approximation  significantly  underestimates  the  probability  by  an  order  of  magnitude.  These  preliminary 
results  provides  the  basis  for  optimism  for  the  proposed  idea. 

The  findings  corresponding  to  this  work  are  disseminated  through  following  publications: 

•  N.  Adurthi,  and  P.  Singla,  “A  Conjugate  Unscented  Transformation  Based  Approach  for  Accurate 
Conjunction  Analysis,”  Richard  Battin  Special  Issue  of  the  AIAA  Journal  of  Guidance,  Control,  and 
Dynamics,  In  Preparation. 

•  N.  Adurthi  and  P.  Singla,  “Principle  of  Maximum  Entropy  for  Probability  Density  Reconstruction: 
An  Application  to  the  Two  Body  Problem,”  2013  Astrodynamics  Specialist  Conference,  Hilton  Head, 
SC,  August  2013. 
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(e)  (x,  z):  upto  3rd  order  (f)  (x,  z):  upto  4th  order 

Figure  5:  MC:  Marginalized  pdf  reconstruction  using  all  3D  moments 
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(d)  (x,  z):  upto  2nd  order 


(c)  (x,  y ):  upto  4th  order 


(e)  ( x ,  z):  upto  3rd  order  (f)  ( x ,  z):  upto  4th  order 

Figure  6:  CUT4:  Marginalized  pdf  reconstruction  using  all  3D  moments 
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(e)  (x,  z):  upto  3rd  order  moments  (f)  (x,  z):  upto  Ath  order  moments 

Figure  7:  CUT6:  Marginalized  pdf  reconstruction  using  all  3D  moments 
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Figure  8:  Snapshots:  Probable  collision  in  the  tail  of  the  density  function. 
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Probability  of  collision  for  threshold  R  Probability  of  collision  for  threshold  R  Probability  of  collision  for  threshold  R  Probability  of  collision  for  threshold  R 


(a) 


(b) 


(c)  (d) 


(g) 


(h) 


Figure  9:  Probability  of  Collision  Computed  Through  Various  Methods 
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4  Efficient  Conjunction  Analysis: 


Interactions  of  and  with  space  debris  and  Resident  Space  Objects  (RSOs)  has  emerged  as  a  significant 
danger  to  most  space  operations.  Current  technology  allows  for  the  tracking  of  debris  down  to  approxi¬ 
mately  five  centimeters  and  from  this,  a  catalog  of  over  21,000  objects  currently  in  orbit  has  been  generated. 
Because  of  the  large  number  of  RSOs  and  the  limited  number  of  sensors  available  to  track  these  objects,  it 
is  impossible  to  maintain  persistent  surveillance  on  all  objects,  and,  therefore  there  is  inherent  uncertainty 
and  latency  in  the  catalog.  Although  small,  the  high  velocities  obtained  by  the  debris  in  low-earth  orbit 
(LEO)  and  geosynchronous  orbit  (GEO)  warrants  a  risk  of  high  damage  to  satellites  and  other  man-made 
objects.  The  examination  of  this  risk,  typically  referred  to  as  conjunction  analysis,  is  a  computationally- 
intensive  process  due  to  the  amount  of  RSOs  considered.  A  recent  incident  in  February  2009  involving  an 
unintentional  collision  between  Russia’s  Cosmos  2251  satellite  and  a  US  Iridium  satellite  clearly  illustrate 
the  need  for  effective  conjunction  analysis.  This  collision  resulted  in  over  500  pieces  of  debris  which  pose 
an  additional  risk  to  other  space  assets.  Brute  force  conjunction  analysis  algorithms  has  combinatorial  com¬ 
putational  complexity  and  hence  require  0(N 2)  computations,  where  N  is  the  number  of  RSOs  considered, 
and  typically  rely  on  numerical  simulations  of  a  mathematical  model  of  the  underlying  physics  of  the  orbital 
debris.  These  simplified  methods,  however,  do  not  account  for  uncertainty  in  the  location  of  the  debris. 

Standard  conjunction  analysis  involves  the  acquisition  of  several  quantities  including  the  time  of  closest 
approach,  the  relative  distance  between  the  objects  in  question,  and  the  probability  of  collision.  A  prelimi¬ 
nary  step  is  often  taken  to  rule  out  any  conjunctions.  In  this  case,  the  geometric  location  of  two  objects  is 
used  to  show  that  the  objects  will  never  cross  paths.  Once  certain  collisions  are  ruled  out,  the  remaining  pos¬ 
sibilities  can  be  investigated.  The  objects  in  question  are  propagated  through  a  dynamical  model,  and  their 
relative  distance  computed.  An  important  step  in  determining  potential  collisions  is  to  delineate  between 
the  orbit  of  each  RSOs.  Objects  are  separated  by  their  orbital  height  from  the  surface  of  the  Earth.  In  this 
sense,  it  is  assumed  that  objects  in  low-Earth  orbit  (LEO)  will  not  collide  with  objects  in  a  geosynchronous 
orbit  (GEO).  The  first  step  is  compute  the  distance  between  the  largest  perigee  value  and  the  smallest  apogee 
value  from  the  two  orbits.  If  this  distance  is  less  than  some  specified  distance  tolerance,  then  the  collision 
of  two  satellites  will  be  examined.  The  next  step  is  to  determine  whether  or  not  a  close  approach  between 
the  two  satellites  occurs.  If  a  close  approach  has  been  determined,  an  error  ellipsoid  must  be  constructed  to 
further  assess  the  risk  of  collision.  The  error  ellipsoid  is  created  about  the  first  satellite  of  interest  using  the 
uncertainty  in  the  satellites  location.  This  results  in  creating  the  ellipsoid  based  on  the  covariance  matrix 
designated  for  the  aforementioned  satellite,  and  ensuring  the  major  axis  is  aligned  with  the  satellite’s  veloc¬ 
ity  vector.  The  next  step  in  analyzing  a  possible  collision  involves  determining  the  probability  of  collision 
of  two  objects. 

Prior  to  computing  the  probability  of  collision,  a  close  approach  between  two  objects  must  be  determined. 
From  this,  the  next  step  is  to  compute  the  relative  velocity  vector  between  the  two  objects  of  interest.  It  is 
assumed  that  all  covariance  data  is  accurate  and  available  at  each  time  step.  For  simplicity,  it  is  also  assumed 
that  the  two  objects  in  question  are  spherical  to  remove  any  attitude  dependencies.  Associated  with  each 
object  is  a  covariance  matrix,  which  is  used  to  specify  an  error  ellipsoid  about  that  object.  Each  error 
ellipsoid  provides  n  standard  deviations  from  the  mean  location  of  the  object,  where  the  value  for  n  is  to  be 
specified  by  the  user  depending  on  the  desired  accuracy.  As  it  is  assumed  that  the  objects  are  independent 
and  identically  distributed  (i.i.d),  the  two  covariance  matrices  are  summed  to  create  a  large  n  —  a  error 
ellipsoid,  where  n  is  the  same  for  both  objects.  This  error  ellipsoid  is  centered  about  the  first  object  of 
interest.  A  virtual  plane  termed  “the  encounter  plane”  is  then  formed  perpendicular  to  the  relative  velocity 
vector,  and  the  error  ellipsoid  is  projected  onto  this  plane,  creating  an  ellipse.  To  represent  all  possible 
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conjunctions,  the  first  and  second  objects  are  merged  to  form  one  spherical  object,  whose  radius  is  equal  to 
the  sum  of  the  two  radii  in  question.  This  combined  object  is  placed  such  that  it  follows  the  relative  velocity 
vector.  From  there,  the  new  object  is  projected  onto  the  encounter  plane.  The  probability  of  collision  is  now 
effectively  reduced  to  a  two-dimensional  computation. 

This  outline  provides  a  complete  brute-force  method  for  determining  conjunctions.  While  complete, 
certain  faults  prevent  this  method  from  its  application  to  all  scenarios.  For  correct  application,  this  method 
requires  that  the  size  of  all  objects  be  known.  Also  required  is  the  condition  that  the  miss  distance  be 
greater  than  or  equal  to  the  combined  object  radius.  The  most  notable  downside  to  this  approach,  however, 
is  the  computational  load  required  to  compute  these  conjunctions.  For  N  objects,  N 2  interactions  must  be 
considered.  While  this  computational  burden  is  the  main  motivation  behind  the  proposed  method,  previous 
attempts  have  been  made  towards  reducing  the  computational  load  of  conjunction  analysis. 

Simple  attempts  such  as  considering  only  objects  that  fall  within  a  specified  orbit,  and  thus  a  certain 
distance  from  of  one  another  have  been  made  in  an  effort  to  reduce  the  computational  load  required  for 
conjunction  analysis  algorithms.  These  efforts,  however,  do  not  take  into  consideration  the  possibility  that 
although  two  objects  may  seem  very  distant  from  one  another,  the  uncertainty  in  their  locations  could  prove 
to  be  high,  warranting  a  possible  collision.  In  this  research  work,  we  introduce  a  hierarchical  approach 
using  a  kd-tree  with  the  aid  of  a  probability -based  distance  metric  known  as  the  Hellinger  distance.  Such 
approaches  have  been  extensively  studied  in  the  gravitational  N-body  problem  literature.  Our  approach  is 
inspired  by  the  great  efficiencies  obtained  by  the  use  of  hierarchical  approaches  and  in  particular  kd-tree 
based  approaches. 

The  nodes  of  a  kd-tree  are  not  physically  created;  they  are  used  simply  to  index  data  points.  The  standard 
kd-tree  implementation  is  described  as  “top-down”  as  this  is  the  direction  of  construction  of  the  tree.  The 
root  box,  which  can  be  described  as  the  space  surrounding  the  entire  data  set,  is  first  initialized.  The  root 
box  is  then  divided  along  a  specified  dimension,  typically  the  longest  available,  creating  two  nodes.  Once 
divided,  the  location  of  each  data  point  is  examined  to  determine  their  respective  placement  in  the  kd- 
tree  -  if  the  particle  is  leftward  of  the  division,  it  is  placed  into  the  leftmost  node,  and  vice  versa  for  a 
rightward  particle.  From  here,  the  procedure  is  repeated  until  only  one  particle  lies  in  each  node  (box). 
Upon  completion  of  the  final  divisions,  the  remaining  boxes,  each  populated  with  one  data  point,  are  termed 
leaves.  This  insertion  procedure  is  termed  recursion,  and  it’s  simplicity  governs  the  favored  implementation 
of  kd-trees.  It  must  be  noted  that  kd-trees  do  not  warrant  the  creation  of  a  physical  mesh  but  rather,  each 
division  is  present  only  to  discern  the  location  of  each  particle.  The  boxes  are  stored  simply  to  index  each 
particle.  It  must  also  be  noted  that  while  the  original  order  of  the  data  set  is  maintained,  during  construction 
of  the  kd-tree  the  data  set  is  organized  in  increasing  order  to  ensure  rapid  construction  of  the  kd-tree.  A 
simplified  kd-tree  algorithm  is  presented  below: 

With  the  kd-tree  implementation  discussed,  the  nearest-neighbor  search  can  be  examined.  A  standard 
nearest-neighbor  (NN)  search  is  conducted  by  computing  the  Euclidean  distance  between  all  possible  com¬ 
binations  of  particles.  By  sorting  the  results,  the  nearest  neighbor  to  each  particle  can  be  determined.  Kd- 
trees  aid  in  NN  searches  by  reducing  the  computational  load  required.  The  computational  complexity  of  a 
standard  “brute-force”  NN  search  is  0(N 2)  whereas  a  kd-tree  based  search  reduces  this  load  to  O(NlogN), 
where  N  is  the  number  of  particles  considered.  A  kd-tree  NN  search  algorithm  is  outlined  below: 

The  standard  NN  search  via  a  kd-tree  implementation  is  clearly  more  efficient  than  a  brute-force  NN 
search  algorithm,  which  requires  the  computation  of  the  distance  between  all  possible  combinations  of 
particles.  It  can  immediately  be  seen  that  the  kd-tree  lends  itself  towards  conjunction  analysis  and  other 
particle  collision  applications. 
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Algorithm  1  Standard  kd-tree  Implementation 
1:  Initialize  size  of  Root  Box  based  on  dataset. 

2:  Reorganize  dataset  based  on  median  value.  Note:  original  order  is  maintained  elsewhere.  Let  p  denote 
a  single  datapoint 

3:  if  p<  division  then  place  p  into  leftmost  node. 

4:  else 

5:  place  p  into  the  rightmost  node. 

6:  if  node  contains  more  than  one  point  then  divide  node. 

7:  else 

8:  store  as  leaf. 

9:  Continue  until  each  node  becomes  a  leaf,  i.e.  each  node  contains  only  one  datapoint. 


Algorithm  2  Kd-Tree  Nearest-Neighbor  Search 
1:  Determine  required  number  of  nearest  neighbors. 

2:  for  p  =  1  to  Number  of  points  do 

3:  begin 

4:  Locate  the  particle  of  interest  in  the  kd-tree. 

5:  Traverse  up  the  tree,  opening  boxes  until  the  required  number  of  nearest  neighbors  is  met. 

6:  Compute  the  distance  from  the  particle  of  interest  to  the  nearest  neighbors  found. 

7:  Store  the  minimum  distance  value,  dmin. 

8:  Compute  the  distance  from  the  particle  of  interest  to  all  remaining  nodes  (boxes)  in  the  kd-tree,  d\)OX. 

9:  if  d\)OX  <  dmin  then  open  the  box  and  compute  the  distance  from  the  particle  of  interest  to  the 

enclosed  particles.  Update  dmin  and  nearest  neighbors. 

10:  end 
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Figure  10:  Comparision  of  Required  CPU  Time  for  Brute-Force  and  Kd-Tree  Nearest  Neighbor  Searches 

It  must  first  be  noted  that  as  each  particle  possesses  an  uncertain  location,  the  mean  value  of  each  particle 
is  used  to  construct  the  kd-tree.  The  main  modification  appears  in  the  NN  search  implementation.  As  op¬ 
posed  to  using  Euclidean  distance,  a  probabilistic  distance  metric  is  to  be  employed  to  compute  the  distance 
between  particles  of  interest.  It  must  be  noted,  however,  that  in  computing  the  distance  between  a  single 
particle  and  a  given  node  (box),  the  Euclidean  distance  metric  is  still  employed.  The  interaction  between 
two  particles,  however,  requires  a  probabilistic  distance  metric.  Here,  the  Kullback-Leibler  divergence,  the 
Bhattacharyya  distance,  and  the  Hellinger  distance  have  been  considered. 

A  preliminary  result  of  this  study  is  depicted  in  Figure  [lOj  showing  the  approximate  CPU  time  for 
both  a  brute-force  and  a  tree-based  nearest  neighbor  computation  as  a  function  of  the  number  of  resident 
space  objects.  It  can  be  seen  that  the  kd-tree  based  approach  offers  a  marked  improvement  in  the  required 
CPU  time  for  the  modified  nearest  neighbor  search  when  compared  to  the  brute-force  approach.  The  kd- 
tree  method  corroborated  results  obtained  by  the  brute-force  method,  ensuring  that  the  Hellinger  distance 
provides  a  link  to  the  actual  probability  of  collision. 

The  findings  corresponding  to  this  work  are  disseminated  through  following  publications: 

•  N.  Adurthi,  and  P.  Singla,  “A  Conjugate  Unscented  Transformation  Based  Approach  for  Accurate 
Conjunction  Analysis,”  2014  AIAA/AAS  Astrodynamics  Specialist  Conference,  San  Diego,  California, 
4  -  7  August  2014. 

•  M.  Mercurio,  “A  Non-Combinatorial  Approach  for  Efficient  Conjunction  Analysis,”  M.S.  Thesis,  Uni¬ 
versity  at  Buffalo,  June  2014. 

•  M.  Mercurio  and  P.  Singla,  “A  Hierarchical  Tree  Code  Based  Approach  for  Efficient  Conjunction 
Analysis,”  2013  Astrodynamics  Specialist  Conference,  Hilton  Head,  SC,  August  2013. 

•  M.  Mercurio,  P.  Singla  and  A.  Patra,  “A  Hierarchical  Tree  Code  Based  Approach  for  Efficient  Con- 
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junction  Analysis,”  2012  AIAA/AAS  Astrodynamics  Specialist  Conference,  Minneapolis,  MN,  Aug. 
2Aug.  5,  2012. 


5  Optimal  Information  Collection: 

Often  the  problem  of  optimally  locating  or  configuring  the  sensor  parameters  is  important  as  it  can 
profoundly  effect  the  performance  of  a  particular  sensor.  For  example  given  the  uncertain  nature  of  space 
debris  and  limited  measurements  or  cost,  one  is  always  interested  in  optimally  configuring/scheduling  the 
sensors  in  future-time  to  make  the  ‘best  measurements’  possible,  as  good  measurements  can  lead  to  better 
decisions  or  improved  estimates.  Often  in  directional  sensors  with  limited  field  of  view  (FOV)  such  as  radar, 
laser,  sonar  or  infrared  ranging  devices,  it  is  important  to  point  the  sensor  at  the  target  (possibly  moving). 
Hence,  the  problem  of  deploying  and  configuring  such  sensors  with  limited  FOV  becomes  an  integral  part 
of  the  measurement  process.  The  concept  of  sensor  management  is  not  new  and  is  still  an  active  field  of 
research  where  many  have  considered  various  utility  functions  to  describe  sensor  performance. 

The  problem  of  optimal  sensor  placement  and  motion  coordination  of  mobile  sensor  networks  has  been 
addressed  in  target  tracking  literature.  The  objective  function  that  measures  the  sensor  performance  is 
generally  taken  as  the  Fisher  Information  matrix  (FIM).  This  sensor  performance  cost  is  apt  as  the  inverse  of 
the  FIM  is  known  as  the  Cramer  Rao  lower  bound  (CRLB)  which  is  the  lower  bound  of  the  covariance  of  the 
parameter  estimates.  As  an  estimate  with  lower  covariance  is  always  desired,  minimizing  this  lower  bound  or 
equivalently  by  maximizing  an  appropriate  norm  of  the  FIM  one  can  achieve  desired  sensor  configurations. 
But  it  is  emphasized  that  the  minimization  is  done  for  a  bound  which  intuitively  implies  better  performance. 
Conventionally,  following  metrics  for  the  fisher  information  matrix  are  used: 

•  D-optimality  (— ln{det(FIM)}) 

•  E-optimality  or  maximum  eigenvalue  (Amax(F/M-1)) 

•  A-optimality  (tr(FIM~1)) 

•  Sensitivity  criterion(— £ r(FIM)) 

Recently,  the  principle  of  maximum  mutual  information  was  used  for  dynamic  sensor  selection  in  case 
of  linear  Gaussian  models.  Often  information  utility  functions  have  neat  analytical  expressions  when  the 
underlying  probability  density  functions(pdf)  are  Gaussian.  For  a  general  pdf,  the  problem  of  computing 
these  information  measures  can  be  a  challenge  and  hence  the  pdf  is  usually  approximated  by  an  equivalent 
Gaussian  pdf.  Gaussian  approximation  in  some  cases  might  not  be  appropriate,  for  example  when  the 
underlying  pdf  is  multi-modal,  a  single  Gaussian  approximation  would  neglect  partly  or  even  completely 
the  high/low  probability  regions. 

The  objective  of  this  research  work  is  to  optimally  manage  various  sensors  configurations  to  tradeoff  be¬ 
tween  competing  indices  such  as  energy  consumption,  coverage  of  domain  and  the  information  gathered.  A 
suitable  information  reward  for  taking  action  (sensor  modality  and  placement)  can  be  the  expected  one  step 
reduction  in  the  total  uncertainty  of  the  system  state.  Due  to  the  uncertain  nature  of  the  target  with  nonlinear 
dynamics,  the  evolution  of  cost  function  is  essentially  stochastic  and  hence  the  problem  can  be  categorized 
as  a  stochastic  optimal  control  problem.  As  the  mutual  information  measure  is  generally  not  convex  in  the 
control  variable,  the  problem  of  even  solving  for  a  numerically  approximate  solution  is  challenging  and  at 
times  intractable. 
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Figure  1 1 :  UAV-Sensor-Target  Schematic 


A  generic  problem  scenario  is  depicted  in  Figure [TT]  where  a  controllable  platform  such  as  an  Unmanned 
Aerial  Vehicle  (UAV)  with  a  configurable  sensor  is  to  move  and  configure  the  sensor  parameters  such  that 
it  can  make  ‘better’  measurements  of  the  uncertain  uncooperative  target.  One  is  also  interested  in  the  cases 
when  there  are  constraints  on  time,  fuel  or  even  the  number  of  possible  measurements.  Let  the  dynamics  of 
the  UAV  with  known  initial  condition  sq  be  described  by 


sk+1  =  F(sk,uk)  (16) 

A  reasonable  assumption  is  made  that  there  is  no  uncertainty  in  the  UAV  state  dynamics.  The  state  of  the 
UAV  typically  consists  of  position,  velocity  and  heading  angle.  The  target  dynamics  with  uncertain  initial 
conditions  xq  ~  p(xq)  is  given  by  an  assumed  non-linear  model 


Xfc+1  =  /(Xfc)  +  vk  (17) 

where  randomness  in  its  motion  is  modelled  by  Gaussian  sequence  vk  ~  J\f{y  :  0,  Q&).  The  measurement 
model  is  a  function  of  the  states  of  the  UAV  (s),  sensor  parameters  ( 6 )  and  the  target  (x)  as  the  measurements 
are  made  with  respect  to  the  coordinate  frame  of  the  UAV  and  are  then  converted  to  the  corresponding 
measurements  in  the  ground  reference  frame. 

zk+ 1  —  h(sk+l,  Xfc+1,  Ok+l)  +  Xfc+1,  (18) 

where  uop+i  ~  :  0,  R/e+i).The  state  pdf  evolves  with  time  according  to  the  Chapman-Kolmogorov 

equation  (CKE).  As  the  measurements  are  made  at  every  time  step  the  pdf  of  the  target  state  is  updated  using 
Bayes’  rule  (BR)  as 


CKE  : 
BR : 


P(Xfc+l  |z&,  Sfc)  — 
p(x/c  +  l  |Zfc+l  7  Sfc) 


/^(xjfe+iIxfe^XjfelzifejSjfe)  dxfc 

_ P{^k-\- 1  l^fc  +  l  )p(-^fc  +  l  | Zfc 

P(zk-\- 1 5Sfc  ) 


(19) 


The  objective  is  to  find  a  sequence  of  control  inputs  [uk,  0k\  to  the  UAV  (and  mounted  sensor)  to  minimize 
a  particular  cost  function.  The  open  loop  optimal  control  problem  can  be  framed  as 
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Nt  —  1 

Min  :  J  =  'ip[sNT]+  ^  < 
k= o 


-  7(xfc,Zfc|sfc,0fc)  H-  Sfc  A.Sk  H-  Ufo  Bilk 


Mutual  Information 


constraint  to  : 


so 


U AV  energy 


s/c+i  =  F(sk,Uk) 

p(x/c+ 1)  =  / p(xfe+i|Xfc)p(Xfc) 
^(0*.)  <  0  and  Cu(uk )  <  0 
P(x0) 


(20) 


(21) 


The  cost  function  consists  of  two  parts:  Mutual  Information,  that  is  to  be  maximized  and  the  energy  (fuel) 
consumed  by  the  UAV  is  to  be  minimized.  Ideally  the  optimal  control  sequence  determined  at  the  initial 
time  is  no  longer  valid  right  after  the  first  measurement  update  and  a  new  optimal  control  problem  needs 
to  be  solved.  This  might  be  computationally  expensive  to  perform  at  every  time  step.  Alternatively  it  is 
logical  to  solve  a  new  optimal  control  problem  only  when  there  is  significant  change  in  the  target  state  pdf 
p(x&|z&,  Sfc)  at  time  k  when  compared  to  the  predicted  target  pdf  p(x&)  at  time  k  for  which  the  optimal 
control  sequence  was  solved. 


DfCL  (p(x/c  |z&5  Sfc)  I  |p(Xfc))  ^  Dfh reShold 


(22) 


In  this  work,  the  principle  of  Dynamic  Programming  (DP)  is  adopted  to  solve  for  the  optimal  action/control: 


J*NT^NT)  =  fp[sNT] 

Jk(sk)  =  min{-J(xfc,zfc|sfc,0fe)  +  ukBuk  +  Jk+1(Fk(sk,uk))} 

uk 


Due  to  the  highly  nonlinear  FOV  constraints,  information  measure  can  change  abruptly  and  optimization 
algorithms  which  involve  gradients  are  usually  not  applicable.  A  discrete  DP  formulation  is  highly  advan¬ 
tageous  as  the  process  does  not  involve  gradients  and  can  also  be  made  parallel.  In  addition  it  is  known 
that  the  computation  time  strongly  depends  only  on  the  level  of  discretization.  Thus  by  adding  multiple 
targets  and  multiple  sensors  the  computation  time  is  weakly  influenced.  Both  Adaptive  Gaussian  Mixture 
Model  (AGMM)  and  Conjugate  Unscented  Transformation  (CUT)  are  used  to  solve  the  CKE  equation  for 
the  optimal  control  problem. 


For  the  purpose  of  simulation,  the  sensor  model  equations  as  described  in  (p~8])  is  composed  of  the  true 
sensor  model  /i(xfc+i,  s^+i,  Ok+i)  and  an  assumed  function  c/(x&+i,  Sfc+i,  9k+i)  that  penalises  the  measure¬ 
ment  made  outside  the  FOV.  Fig[l2fc,  shows  one  kind  of  penalty  g(.). 

In  this  case  study,  the  optimal  trajectory  to  track  a  moving  target  is  sought.  For  this  simulation  the  target 
is  assumed  to  be  an  airplane  with  state  vector  x  =  [£,  77, 77,  making  a  Coordinate  turn  (CT): 


'1  sin(QT) 

1  Q 

0  cos(QT) 

PI  l—cos(QT) 
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%k- 1  +  vk- 1, 


(23) 


[^,  rj\  are  the  position  coordinates  of  the  airplane  and  [£,  77]  are  the  corresponding  velocities,  Q  is  the  turn 
rate.  Here  the  T  =  Is  is  the  discrete  time  step  for  the  target  dynamics.  The  assumed  process  noise  Vk-\ 
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Figure  12:  Sensor-Target  Schematic 


has  covariance  matrix  Qk-i  with  parameter  (Li,  L^)  =  (0.02, 10  5).  The  measurement  model  consists  of 
range  and  bearing  angle  measurements  within  the  FOV  of  the  sensor. 


/i(x,  s,  6)  = 


v(f~  Si)2  +  (v  -  S2)2 


tan 


(24) 


The  FOV  consists  of  a  circular  sector  with  angle  60°  and  radius  50m.  A  measurement  update  is  performed 
only  when  the  true  target  position  falls  within  the  FOV.  The  constant  measurement  noise  covariance  is 
Rk+i  =  diag[ 0.52  m2,  (l7r/180)2  rad2].  As  the  measurement  model  is  not  a  function  of  the  velocity 
of  the  target,  the  mutual  information  can  be  computed  from  the  marginalised  pdf  .The  DP  algorithm  is 
initially  computed  for  the  whole  simulation  time  i.e  25  time  steps.  A  new  optimal  trajectory  is  computed 
at  time  step  12  when  the  threshold  D threshold  in  eqn  ([22])  exceeds  10  by  solving  the  DP  problem  from 


time  step  4  =  12  to  trf  =  25.  To  calculate  the  information  and  propagate  the  state  pdf  the  minimal 
cubature  points  of  CUT6  are  used.  For  this  particular  simulation  the  initial  condition  uncertainty  is  taken 
as  p(x 0)  =  0.57V(x  :  /ii,  Pi)  +  0.5A/’(x  :  fi 2,  P2),  where  the  Gaussian  components  represent  information 
from  two  different  sources.  The  values  for  the  means  and  covariance  are: 


fi±  =  [10,  -1, 10,  5,  -5tt/180]t  /ii  =  [10,  2, 10, 3,  -3tt/180]t 
Pi  =  diag[  1, 0.5, 1, 0.5,  lO”3]  P2  =  diag[  1, 0.5, 1, 0.5,  lO”3] 


A  random  point  from  this  distribution  is  selected  as  the  true  initial  position  of  the  target.  It  is  emphasized 
that  this  true  location  along  with  its  trajectory  is  unknown  to  the  DP  problem  and  is  only  used  to  generate 
random  measurements  when  the  target  is  in  the  sensor  FOV,  otherwise  there  is  no  measurement  update  of 
the  state  pdf.  Fig  |13(f)|  shows  the  drop  in  the  joint  entropy  of  the  target  state  pdf  with  every  measurement 
compared  to  the  pure  propagation  of  the  target  pdf.  During  the  initial  time  steps  it  can  be  seen  that  there  is 
no  decrease  in  the  uncertainty  as  the  target  pdf  evolves.  This  is  because  the  target  is  not  visible  to  the  sensor. 
Figs  1 1 3 (a)]  - 1 1 3 (e) |  show  snapshots  of  the  simultaneous  motion  of  the  sensor  along  the  optimal  trajectory  and 
the  updated  state  pdf  of  the  target. 


Finally,  we  solve  an  information  optimization  problem  to  manage  a  network  of  ground  based  sensors 
to  maximize  the  mutual  information  on  the  system  of  space  objects  being  observed  by  the  network.  A 
limited  look-ahead  policy  is  used  to  solve  the  dynamic  programming  problem  that  schedules  the  elements  of 
a  ground  based  surveillance  network  to  maximize  the  mutual  information.  Simulation  results  will  be  used 
to  validate  the  sensor  management  algorithms  reported  in  the  paper. 

Figure  |T4]  shows  various  time  snapshots  of  the  sensor  management  algorithm  performance  where  10 
RSOs  are  tracked  using  5  radar  observation  stations  with  known  sensor  uncertainties  at  various  locations  on 
the  Earth.  Unscented  transformation  is  used  to  evaluate  the  FIM  and  computation  of  the  mutual  information 
performance  measure  for  sensor  management.  Figure  [15]  shows  the  time  history  of  the  Frobenius  norm  of 


22 


x-axis 

(a)  Time  step  k=4 


(c)  Time  step  k=13 


0  20  40  60  80  100 

x-axis 


(e)  Time  step  k=22 


100 

80 

»  60 


0  20  40  60  80  100 

x-axis 


0  20  40  60  80  100 

x-axis 


(d)  Time  step  k=18 


(f)  Decrease  in  Joint  Entropy 


Figure  13:  Optimal  Trajectory 


the  covariance  matrix  of  each  object  under  surveillance  using  the  optimal  sensor  schedule  generated  by  our 
approach.  The  reduction  of  covariance  of  all  the  objects  being  tracked  demonstrates  the  optimism  of  the 
technical  approach  discussed  in  the  full  paper.  The  full  paper  presents  all  the  technical  details  of  the  sensor 
management  algorithms. 

The  findings  corresponding  to  this  work  are  disseminated  through  following  publications: 

•  N.  Adurthi,  P.  Singla  and  M.  Majji,  “Optimal  Sensor  Tasking  for  Space  Situational  Awareness,” 
Richard  Battin  Special  Issue  of  the  AIAA  Journal  of  Guidance,  Control,  and  Dynamics,  In  Prepa¬ 
ration. 

•  N.  Adurthi,  P.  Singla  and  M.  Majji,  “Information  theoretic  Optimal  Sensor  management  for  Efficient 
Space  surveillance,”  2014  AIAA/AAS  Astrodynamics  Specialist  Conference,  San  Diego,  California,  4 
-  7  August  2014. 

•  N.  Adurthi  and  P.  Singla,  “Information  Driven  Optimal  Sensor  Control  for  Efficient  Target  Localiza¬ 
tion  and  Tracking,”  2014  American  Control  Conference,  Portland,  OR,  June  4-6,  2014. 

•  N.  Adurthi,  P.  Singla  and  T.  Singh,  “Optimal  Information  Collection  for  Nonlinear  systems-  An  Appli- 
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Figure  14:  Snapshots  of  the  simulation  at  various  times 
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Figure  15:  Reduction  in  Frobenius  norm  of  covariance  of  each  object  with  time 

cation  to  Multiple  Target  Tracking  and  Localization,”  2013  American  Control  Conference,  Washington 
D.C. 

•  R.  Madankan,  R  Singla  and  T.  Singh,  “Optimal  Information  Collection  for  Source  Parameter  Estima¬ 
tion  of  Atmospheric  Release  Phenomenon,”  2014  American  Control  Conference,  Portland,  OR,  June 
4-6,  2014. 

•  K.  Vishwajeet,  N.  Adurthi  and  P.  Singla,  “Optimal  Information  Collection  for  Space  Situational 
Awareness,”  2014  SIAM  Conference  on  Uncertainty  Quantification,  Savannah,  GA,  March  31 -April 
3,  2014. 

•  N.  Adurthi,  R.  Madankan  and  P.  Singla,  “Optimal  Information  Trajectory  Design  for  Dynamic  State 
Estimation,”  2014  SIAM  Conference  on  Uncertainty  Quantification,  Savannah,  GA,  March  31-April 
3,  2014. 

6  Polynomial  Chaos  based  Bayes  Filter: 

Two  new  recursive  approaches  have  been  developed  to  provide  accurate  estimates  for  posterior  moments 
of  both  parameters  and  system  states  while  making  use  of  the  generalized  Polynomial  Chaos  (gPC)  frame¬ 
work  for  uncertainty  propagation.  The  main  idea  of  the  gPC  method  is  to  expand  random  state  and  input  pa¬ 
rameter  variables  involved  in  a  stochastic  differential/difference  equation  in  a  polynomial  expansion.  These 
polynomials  are  associated  with  the  prior  pdf  for  the  input  parameters.  Later,  Galerkin  projection  is  used 
to  obtain  a  deterministic  system  of  equations  for  the  expansion  coefficients.  The  first  proposed  approach 
(gPC-Bayes)  provides  means  to  update  prior  expansion  coefficients  by  constraining  the  polynomial  chaos 
expansion  to  satisfy  a  specified  number  of  posterior  moment  constraints  derived  from  the  Bayes’  rule.  The 
second  proposed  approach  makes  use  of  the  minimum  variance  formulation  to  update  gPC  coefficients.  The 
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main  advantage  of  proposed  methods  is  that  they  not  only  provide  point  estimate  for  the  state  and  parameters 
but  they  also  provide  statistical  confidence  bounds  associated  with  these  estimates. 

To  illustrate  the  effectiveness  of  the  developed  ideas,  let  us  consider  the  example  involving  the  the  Duffing 
oscillator: 


(25) 

(26) 

For  simulation  purposes,  nominal  parameter  values  are  assumed  to  be  given  as: 

r]  =  1.3663,  a  =  -1.3761  /3  =  2 


x  +  r]x  +  ax  +  /3x 3  =  sin(3£) 
x(tk) 


y  (*k)  = 


x{tk) 


+  vk 


The  initial  states  are  assumed  to  be  normally  distributed: 

x(0)  =  J\f(x0 1  -  1,  0.25),  £(0)  =  M{x 0|  -  1, 0.25) 

Hence,  'ipk(£)f s  are  chosen  to  be  Hermite  polynomials  to  describe  Gaussian  distribution  of  states.  As  well, 
4th  order  gPC  expansion  is  considered  to  analyze  the  effect  of  initial  condition  uncertainty. 

To  corroborate  the  efficacy  of  the  PCQ  approach  to  capture  the  evolution  of  the  statistics  of  states  of  ([25]), 
relative  error  in  Frobenius  norm  of  the  difference  between  different  moments  of  states  with  respect  to  105 
Monte  Carlo  runs  at  t  —  2sec.  is  evaluated.  Table [3]  shows  that  the  relative  error  decreases  as  the  number  of 
quadrature  points  increases.  It  is  clear  that  one  can  obtain  a  better  approximation  for  three  central  moments 
using  only  16  quadrature  points,  relative  to  the  103  MC  runs. 


Table  3:  Relative  error  in  Frobenius  norm  of  the  difference  between  moments  of  states  and  105  Monte  Carlo 
runs  at  t  =  2sec. _ 


Number  of  Quadrature  Points 

Mean 

2nd  Central  Moment 

3rd  Central  Moment 

l2 

4.5526% 

100% 

100% 

22 

0.3217% 

20.3050% 

98.3149% 

32 

0.0329% 

3.9559% 

28.5219% 

42 

0.0537% 

0.5202% 

2.5084% 

103  MC  Simulations 

0.1199% 

6.0715% 

99.2219% 

To  verify  the  efficiency  of  our  method,  we  compared  the  performance  of  the  proposed  methods  with 
the  extended  Kalman  Filter  (EKF)  and  Particle  Filter  (PF)  results.  The  measurement  data  is  assumed  to 
be  available  at  a  sampling  frequency  of  1Hz.  A  random  sample  of  initial  conditions  is  taken  from  initial 
condition  distribution  to  generate  the  noise-free  measurement  data.  The  noise-free  measurement  data  is  then 
corrupted  with  a  Gaussian  white  noise  of  zero  mean  and  variance  being: 


R  = 


a2  0  \ 
0  a2  ) 


a  is  assumed  to  be  0.05  in  our  simulations. 

Fig.  1 16(a)]  and  Fig.  1 16(b)]  illustrate  the  state  estimation  error  for  x  and  x  by  using  EKF  method,  respec¬ 
tively.  The  solid  blue  line  represents  the  difference  between  the  true  value  and  its  mean  estimate.  Dashed 


26 


(a)  x{t)  (b)  x(t) 

Figure  16:  Error  and  3cr  Bounds  for  the  EKF  Approximated  Posterior  Mean  for  Duffing  oscillator 


(a)  x(t)  (b)  x(t) 

Figure  17:  Error  and  3a  Bounds  for  PF  Approximated  Posterior  Mean  for  for  Duffing  oscillator 


green  line  shows  —3a  bound  while  the  dashed  red  line  represents  the  3a  bound.  From  these  plots,  it  is 
clear  that  the  state  estimation  error  increases  significantly  with  time  although  it  is  always  bounded  by  3a 
bounds.  The  poor  performance  the  EKF  can  be  attributed  to  strong  nonlinearities  and  sparse  data  resulting 
from  sampling  at  1  Hz. 


The  state  estimation  error  for  x  and  x  by  using  Particle  Filter  have  been  shown  Fig.  1 17(a)] and  Fig.  1 17(b) 


respectively.  The  solid  blue  line  represents  the  difference  between  the  true  value  and  its  mean  estimate. 
Dashed  green  line  shows  minimum  bound  while  the  dashed  red  line  represents  the  maximum  bound.  These 
plots  show  that  the  state  estimation  error  decreases  during  the  time,  while  using  PF. 


Furthermore,  Fig.  [18]  shows  the  error  in  state  estimates  along  with  its  3a  bounds  using  the  gPC  based 
minimum  variance  estimator.  Once  again,  the  estimation  error  along  with  3a  bounds  converge  to  zero  over 
the  time  which  can  be  again  attributed  to  the  posterior  density  function  being  a  delta  function  as  number  of 
measurements  increases. 


Fig.  [19]  shows  the  error  in  state  estimates  using  the  gPC-Bayes  method  for  various  values  of  Nm.  The 
solid  blue  line  represents  the  difference  between  the  true  value  and  its  mean  estimate.  Dashed  green  line 


27 


(a)  x(t) 


(b)  x(t) 


Figure  18:  Error  and  3 a  Bounds  for  the  Minimum  Variance  Approximated  Posterior  Mean  for  Duffing 
oscillator 


Table  4:  RMSE  error  in  mean  estimate  of  states  x  and  x  while  using  different  estimation  methods 


Method 

&X 

EKF 

0.1307 

0.3858 

PF 

0.0408 

0.0508 

minimum  variance 

0.0359 

0.0527 

Nm  =  1 

0.1173 

0.0695 

gPC-Bayes 

Nm  =  2 

0.0347 

0.0531 

Nm  =  3 

0.0336 

0.0528 

and  dashed  red  line  represent  the  min  and  max  bounds  on  estimation  errors,  respectively.  It  is  clear  that 
estimation  error  and  corresponding  3cr  bounds  for  estimation  error  converge  to  zero  over  the  time.  This  is 
due  to  the  fact  that  posterior  density  function  finally  converges  to  a  dirac-delta  function  around  the  truth 
which  is  expected  as  number  of  measurements  increases  over  the  time.  Also,  it  should  be  noticed  that  3a 
bounds  becomes  more  and  more  tighter  as  one  increases  the  number  of  matching  moment  constraints,  i.e., 
Nm.  From  these  results,  it  is  clear  that  the  proposed  methods  perform  very  well  in  not  only  estimating  the 
posterior  mean  but  posterior  density  function  also. 

To  summarize,  the  root  mean  square  error  over  time  between  the  mean  estimate  of  states  and  their  true 
value  has  been  shown  in  Table  [4]  As  this  table  represents,  gPC-Bayes  and  PF  method  perform  very  well 
in  estimation  of  both  states  x  and  x,  while  EKF  results  in  high  error  between  the  mean  estimate  and  actual 
value  of  the  states.  It  is  clear  from  Table  [4]  that  by  increasing  the  number  of  matching  moment  constraints 
(Nm)  in  gPC-Bayes  method,  the  error  in  estimation  of  states  decreases. 

The  findings  corresponding  to  this  work  are  disseminated  through  following  publications: 

•  R.  Madankan,  P.  Singla,  T.  Singh  and  P.  Scott,  “Polynomial  Chaos  Based  Method  for  State  and  Pa¬ 
rameter  Estimation,”  AIAA  Journal  of  Guidance,  Control  and  Dynamics,  Vol.  36,  No.  4  (2013),  pp. 
1058-1074,  July  2013,  DOI:  10.2514/1.58377, 

•  R.  Madankan,  P.  Singla,  T.  Singh  and  P.  Scott,  “Polynomial  Chaos  Based  Method  for  State  and  Pa¬ 
rameter  Estimation,”  2012  American  Control  Conference,  Montreal,  Canada,  June  27-29,  2012. 
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(a)  Estimation  Error  for  x  (Nm  =  1) 


(b)  Estimation  Error  for  x  ( Nm  =  1) 


(c)  Estimation  Error  for  x  (Nm  =  2)  (d)  Estimation  Error  for  x  (iVm  =  2) 


(e)  Estimation  Error  for  x  (Nm  =  3)  (f)  Estimation  Error  for  x  (Nm  —  3) 


Figure  19:  Estimation  Error  and  3cr  Bounds  for  the  gPC-Bayes  Approximated  Posterior  Mean  for  for  Duffing 
oscillator 
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7  Concluding  Remarks 


In  summary,  significant  progress  has  been  made  towards  characterizing  non-Gaussian  state  density  func¬ 
tion.  Two  new  methods  adaptive  Gaussian  mixture  model  (AGMM)  and  conjugate  unscented  transformation 
(CUT)  have  been  developed  for  this  purpose.  AGMM  method  solves  the  Fokker-Planck-Kolmogorov  equa¬ 
tion  associated  with  orbital  dynamics  model.  Furthermore,  sparse  approximation  tools  have  been  used  to 
identify  the  Gaussian  kernels  which  best  approximate  the  state  pdf.  The  CUT  methodology  provides  effi¬ 
cient  means  to  compute  higher  order  moments  of  state  density  functions.  The  primary  objective  of  CUT 
methodology  is  to  find  a  fully  symmetric  sigma/cubature  point  set  with  reduced  number  of  points  that  is 
equivalent  to  the  set  of  cubature  points  of  Gaussian  quadrature  product  rule  of  same  order.  Equivalent  to 
same  order  implies  that  for  a  polynomial  of  order  2 m  —  1  in  generic  7V-dimensions,  both  the  new  reduced 
sigma  point  set  from  the  proposed  method  known  as  Conjugate  Unscented  Transform  method  (CUT)  and  the 
mN  quadrature  points  from  the  Gaussian  quadrature  product  rule  result  in  same  order  of  relative  percentage 
error.  In  this  work,  a  closed  form  expression  for  these  new  sets  of  point  is  provided  to  satisfy  up  to  8  central 
moments.  It  is  shown  that  the  proposed  method  provides  a  significant  reduction  in  function  evaluations  to 
compute  multi-dimension  expectation  integrals.  For  example,  the  proposed  method  needs  only  355  and 
745  function  evaluations  to  compute  the  expectation  integral  for  a  polynomial  function  of  degree  8  in  the 
5-  and  6-dimensional  space,  respectively  whereas  the  Gaussian  quadrature  product  rule  would  need  3,125 
and  15,625  function  evaluations  for  the  same  5-  and  6-dimensional  space  respectively.  Finally,  optimal 
control  problem  is  posed  to  optimize  sensor  locations  and  other  modalities  to  better  track  a  target  object. 
The  main  highlight  of  this  work  is  to  compute  information  theoretic  metrics  corresponding  to  non-Gaussian 
target  state  pdfs  to  describe  the  current  situation  of  the  target  state. 
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Training  (Be  Specific) 

Supplies 

Other  Expenses  (Be  Specific) 

Total  Resource  Requirements 

Report  Document 
Appendix  Documents 
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Response  Location 


Country: 

United  States 

Region: 

NY 

City: 

Buffalo 

Postal  Code: 

14260 

Long  &  Lat: 

Lat:  42.768398284912,  Long:-78.887100219727 

